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I provide a transfer matrix method for the Foldy-Wouthuysen representation of the Dirac equation. 
I develop the relationship between the reflection and transmission coefficients of the Dirac spinors 
and the wavefunction in the transformed representation. I develop a WKB approximation for Dirac 
fermions which has the same elegant form as the WKB solution to Schrédinger’s equation. My WKB 
approximation is to all orders and includes the semi-classical turning point. I provide an extension to 
fully 2-dimensional periodic structures by Fourier methods for band engineering. I verify my methods 
for all energies by comparison with analytic solutions developed in the Dirac spinor representation. I 


also include a rich appendix detailing my research into the Green’s functions of Dirac fermions. 


1. Introduction 


In Refs. [1, 2, 3] the pseudo-relativistic dispersion of 
2-dimensional (2D) materials was investigated theoretically 
for the emergence of extra Dirac points with 1D periodic bar- 
riers. The method employed was the application of Bloch’s 
theorem to transfer matrices of Dirac spinors. These band 
engineering studies attracted considerable interest since the 
vanishing band gap of 2D materials is an Achilles heel for ap- 
plications. Other methods for band engineering that open the 
way to creating 2D semiconductor materials include strain 
engineering [4], vertical stacking [5], and chemical doping 
[6]. 

Armed with the knowledge to create semiconductor 2D 
materials a theoretical study of the operational window for 
negative differential resistance was pursued by other authors 
[7]. The mathematical method employed in their study was 
the application of transfer matrices to the Dirac spinors of a 
rotated Dirac Hamiltonian derived in Ref. [8]. 

Motivated by the study of unconventional transmission 
and density of states, the modeling of Dirac fermions has 
moved beyond unrealistic square barrier systems. In exper- 
imental devices barriers are more often smoothly varying, 
semi-classical (WKB) methods are needed. An approach us- 
ing separate Hamiltonians for electrons and holes in Ref. [9] 
gave way in Ref. [10] to an expansion in powers of hf partly 
resembling the WKB method for the Schrédinger equation. 
However it is acknowledged by the authors of Refs. [9, 10] 
that their mathematical approaches to semi-classical analy- 
sis can be divergent at the classical turning point. 

The contribution to the theory of tunneling in 2D ma- 
terials that I present in this article is the application of the 
Foldy-Wouthuysen (FW) representation of the Dirac equa- 
tion. 

Whilst the FW transformation was first conceived as a 
semi-relativistic representation of the Dirac equation [11], 
I will show, through my rigorous derivations in Section 2 


“Corresponding author 


& doostmb@gmail.com, https://gofund.me/34294631 (M.B. Doost) 
ORCID(S): 0000-0003-4682-2889 (M.B. Doost) 


to Section 13 and my analytic and numerical verification of 
Section 14 and Appendix B, that the FW representation re- 
mains an exactly accurate representation of the Dirac equa- 
tion for all energies. 

I find that the FW representation allows a derivation of 
the relativistic WKB approximation, Section 7 to Section 10, 
in the same elegant form as the approximate WKB solution 
to Schrédinger’s equation. Since at low energies the FW 
representation reduces to the Schrédinger equation, I have 
found, in Section 11, a simple method to mathematically de- 
scribe the connections between WKB regions. 

Fourier analysis is used in Section 12 to apply my ap- 
proach to fully 2D periodic structures. Section 12 offers 
the prospect of extending the scope of the band engineering 
studies made in Refs. [1, 2, 3]. My numerical and analytic 
results showing that the FW equation is an exactly accurate 
representation of relativistic fermions, for all energies, will 
be important here. The Fourier transform method in Section 
12 requires expansion of the wavefunction in plane waves of 
momentum beyond semi-relativistic energies. 

In Appendix E I make the reader aware of a limit to 
appealing for solutions of the Dirac equation from its FW 
representation, I reveal and discuss a Green’s function (GF) 
paradox arising when transforming between the two repre- 
sentations. 

My article is organized as follows: Section 2 outlines the 
derivation of the FW representation from the Dirac equa- 
tion. Section 3 gives the transmission of the Dirac equation 
in terms of the transmission of its FW representation. Sec- 
tion 4 derives the boundary conditions of the FW equation 
at sharp steps in potential. Section 5 and Section 6 develop 
the transfer matrices for square barriers and delta potentials. 
Section 7, Section 8, Section 9 and Section 10 give a 1st, it- 
erative and 2nd order WKB approximations for two interpre- 
tations of the FW equation. Section 11 gives the connecting 
formulae between regions where the WKB approximation 
is appropriate. Section 12 extends my approach to include 
periodic structures. Section 13 details how I evaluated the 
boundary conditions. Section 14 gives numerical results for 
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examples of tunneling through barriers and resonant diodes. 

The appendices are organized as follows: Appendix A 
analytically justifies my commutation of FW wavefunction 
operators. Appendix B analytically calculates the reflection 
at a step in both the Dirac spinor and FW representations for 
comparison. Appendix C examines the tunneling of a mass- 
less fermion through a magnetic delta barrier in the FW rep- 
resentation. Appendix D calculates GFs of the FW equation. 
Appendix E uncovers a relativistic GF paradox. 


2. Dirac equation and its Foldy-Wouthuysen 
representation 


The Dirac equation for a scalar potential V is given in 
Hamiltonian form by 


[ca - p+ Bmc?| ¥ = HY =(E+V)¥, (1) 


a +1 0 
= (2, aR iA a) m 


with Oy yz = VN xy,2 and p= y,. 


a) 2| (5) 


p = -ihV = -ih | —, —, — 
4 ' ; E Oy Oz 


is the momentum operator for the electron which is of mass 
m and energy E. The charge and current density are given 
by the following well known expressions: 


p=Y, (6) 


jo Pov. (7) 


There are two linearly independent solutions of the free par- 
ticle Dirac equation, 


Fmc* — E-V 
(4) _ 0 IP,X + Ipyy 
zy = 0 exp ( h ? (8) 
c (py + ipy) 
0 
=mc2 —- E— ip,x +ipyy 
wee) Fm EE Te “ — ). (9) 


¢ (Dy + ipy) 
0) 


each can be normalised to one electron per unit volume by 
Eq. (6). 

By making use of the FW transformation I will decouple 
the differential equations of the Dirac equation to give four 
independent identical equations. I will show that by calcu- 
lating the propagation for only one of these components I 
may obtain the full GF and transmission. 

Under unitary transformation the Dirac equation becomes 


et!S eS oti Swe (E+ V) et Sh, (10) 
Foldy and Wouthuysen [11] found that their rotation 
tS = cos + fa sind, (11) 
Pp 


transformed the Dirac Hamiltonian into diagonal and anti- 
diagonal parts 


B (mc? cos (20) + ¢ |p| sin (26) 


+ca-p (cos (20) — a sin (20) 
p 
=et!SHe-iS (12) 


Subsequently it was chosen [11] 


2 
cos (26) = ——“{___, (13) 
Ve2p2 + m2c4 
sin (26) = CA (14) 


fc2p2 + mec4 
so that by direct substitution of Eq. (13) and Eq. (14) into 
Eq. (12) 


[ove + met -1V +E] F=0. (15) 


Eq. (15) shows there are two possibilities for the FW equa- 


tion 
EVP em-V+D)F=0, 16) 


yo) corresponds to the electron, Y) to the positron. The 
rotated wavefunction P is related to ¥ by 


Joos + pa sino] v=, (17) 
P 
Taylor expanding +/ p2c? + m2c* in Eq. (16) gives 

[FL - (V +E +mc’)| PO =0, (18) 
where 


pene 1[2) .2(2) 42127. (19) 
2 \)mc 8 | mc 16 | mc aie Ue 


Eq. (16) and Eq. (18) are two interpretations of the FW equa- 
tion. 

I demonstrate the unsuitability of my FW approach for 
calculating magnetic barrier tunneling in Appendix C. I dis- 
cuss the GFs of the FW Eq. (18) in Appendix D. I uncover 
and discuss a GF paradox of the FW Eq. (18) in Appendix E. 
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3. Relationship between the transmission of 
the Dirac equation and its 
Foldy-Wouthuysen representation 
Throughout my numerical and analytic results, Section 


14 and Appendix B, I will calculate the reflection coefficient 
R to deduce the transmission coefficient T as 


T=1-R. (20) 


Calculating R, rather than T directly, does not require the 
normalisation of the wavefunction in both the incident and 
transmitted regions. 

Let ¥ be given by Eq. (8) and Eq. (9) as 


ip,x +ipyy 
W = uexp (A) (21) 
h 
taking of the form 
eos ip,x +ipyy 
W = trexp (7), (22) 


I find by application of Eq. (17) to Eq. (21) 
i, =u, cos0 + (d, —id,) using, 
Uy = u,cos6 + (d, +id,) u; sind, 
is 2 ( x 5) 3 (23) 
U3 = U3 COS 8 — (d,. - id,) uy sind, 
ty = uycos0 — (d, +id,)u, sind, 

where d = p/|Al. 

Eq. (21), Eq. (22), and Eq. (23) show ¥ in exponential 
planar form which is useful for calculating the fermion prop- 
agation through homogeneous space. 

The incident (i) and reflected (7) waves for planar struc- 
tures in the Dirac equation description are related to the re- 


flection coefficient in the well known way: 
Bry 
goes. (24) 
(P. pr 


For the two linearly independent solutions of the Dirac equa- 
tion, Eq. (8) and Eq. (9), Eq. (24) gives rise to two degenerate 
expressions for the reflection coefficient: 


(r) 
(MiP + lal?) 
R, = ————____, (25) 


i) 
(Mi? + Nal’) 


(r) 
(2)? + ['¥517) 
————————— er (26) 


(i) 
(Ia)? + [¥51") 


I take note of Eqs. (23) linking w to u and calculate when d 
is real 


lz? oF |ii,|7 = [u;|? + |u|” , (27) 


lita]? + |ais|” = [ue |? + |usl’ (28) 


therefore R described by the Dirac equation can be rewritten 
in terms of the FW wavefunctions: 


(Zee 


R= (29) 


Gal r ef) ‘ 
(ieee 


“Pa ° 
2 3 


Since the components of ® are degenerate I am only re- 
quired to evaluate the propagation of a fermion by the FW 
Eq. (18) once when calculating R and T. 


4. Boundary Conditions of the 
Foldy-Wouthuysen equation at a sharp step 
Consider a sharp step in potential at x = a. I will inte- 


grate the FW Eq. (18) through the boundary, normal to the 
boundary 


; a —~ 41) AX 
1 FL-(V+E *\| pH =" =0. Bl 
ima | [= ( + mc )| = 0. (31) 


In this section I introduce the operator £ which com- 
mutes with cp, and is defined by 


LY = —Lep,¥ = cp, LP. (32) 


I give a discussion of Eq. (32) in Appendix A. 
Since 


a4 ~ 
lim / (V+ E+ me?) Go) & =0, (33) 
(a,-a_)>0 Sa ch 


due to the function integrated being finite but the integration 
range tending to zero, I have 


=, a 
[ize| Fe ite / cel® <0. (34) 
a_ (a,-a_)>0 Jq ch 


I can understand Eq. (34) when I consider that 


LY = ae i) 


(35) 


Eq. (34) shows that ¥ is continuous with respect to the op- 
eration of L. 7 

In order to evaluate £'Y for the mathematical description 
of the boundary conditions note 


+Lep,P* = (V+ E+ mc?) B® = scp, LP), (36) 
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therefore 


V+E+mc 


-~ 2~ ee 
LPH = + = yO — sqbyp) | (37) 
x 


For my second boundary condition I note that pe always 
gives a finite observed value, therefore ¥ is continuous ev- 
erywhere. 

My boundary conditions are in stark contrast to the con- 
tinuity conditions of the Klein-Gordon equation, in that case 
continuity is with respect to momentum and wavefunction. I 
verify my boundary conditions numerically and analytically 
in Section 14 and Appendix B. 


5. Transfer matrices for a square potential 
barrier 


Inside a layer thickness D the transfer matrix of W is 
given by 


iDp,. 
i ers 0 ay dg 
iDp ee cal oe 
0 exp (+ Ps ) b, bo 


where dy and a, denote forward traveling FW wave parts 
and Bo and Dy denote backward travelling FW wave parts, 0 
denotes on the left, 1 denotes on the right. The fermion is 
refracted at the angle @ to the layer boundary normal, where 


» (38) 


Px = pcos. (39) 


For a boundary at a sharp step, applying the wavefunc- 
tion continuity condition set out in Section 4 gives 


do + Bo = ay +b, ‘. (40) 


Secondly, applying wavefunction continuity under the oper- 
ation of £, derived in Section 4, to Eq. (40) gives 


ad = agby = aa) am aby f (41) 


Eq. (40) and Eq. (41) can be written in matrix form as 


igs 1-2)\Vle) {- 
: 0 Pra | ee a ee (42) 


a ay ay ~| [x 
(1-2) 141) by bo 
a a 


In the event of tunneling, wave modes become an expo- 
nentially decaying and growing set. Transfer matrix multi- 
plication combines these large and small values, this process 
is numerically unstable and inaccurate. Instead I suggest 
using the method of Refs. [12, 13] to combine my transfer 
matrices. The method of Refs. [12, 13] addresses these in- 
stability issues by separating the exponentially growing and 
decaying terms. 


6. Transfer matrix for a delta potential 
barrier 
Integrating Eq. (18), with V(x) = g6(x) and with respect 
to dx/ch, across [O_, 0, gives 


[Fiz (a +5)" = 3 (a +5) . (43) 
OL ch x=0 
Continuity of the wavefunction gives 
dy + by = G +b,. (44) 


Eq. (43) and Eq. (44) can be written in the following matrix 


form 
i- ig _ ig - 
2cha,, 2cha, a) 
ig ig ~ , 
1 
(+ 2cha, ) ( - 2cha, ) y 


(45) 


iS) 
Sl 


rl 


ay 
o 
So 


where I define G, @,, bg, and 5, in the same way as Section 
3. 


7. 1st order relativistic WKB approximation 


Consider the following result from Appendix D derived 
from the correct fermion interpretation of the FW equation 
in terms of £, 


2 
2 
T(h= ee i) (46) 


2 Im* 
(a? — a?) sin? ( a) + (2a,a,)” 


for a rectangular barrier given by 


_ Jj V~ for |x| <a, 
v= { 0 for |x|>a. om) 


Eq. (46) can be approximated as 


2 
2a; 4im*qa 
raos( ; s) exp (- — ) (48) 
a, — a 


by assuming exponential wavefunction decay within the bar- 
rier. Taking logarithms of Eq. (48), the non-logarithmic term 
dominates the logarithmic term so that 


4im* 
log, T (k) x -——“ 


(49) 


Assuming independence of transmission events, the total 
transmission probability for crossing an overall barrier is 


T;, (50) 


| 
Q 


2im; qi;Ax 


= 2i 
log. [71 * -—— = -5 | pdx, (51) 


i=1 
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where the T; are the transmission probabilities of the indi- 
vidual barriers. I see from Eq. (51) that inside a smoothly 
varying layer the FW wavefunction is approximately given 


by 
N 
és ; 2aj.a 
P ~ ex st fp ax) eV 86152) 
E ( h ' II a2 — a2 


k q 


8. Iterative relativistic WKB approximation 


Consider a wavefunction of the form 


Y= exp (2), (53) 


then the boson FW interpretation Eq. (16) becomes 


2 2 ~ 
fe) — in + met — (V+ E) P=0. 
x 


Ox 
(54) 


Eq. (54) can be exactly rewritten as 


2 2 
(2) - ine +mct=(V+EY. (55) 


I now make use of (E + V)* = m?c* + pc? in Eq. (55) and 
find that as expected 


ty _ 20 ee 
BLY = pL, oe _ reece (56) 
I know from Section 7 
o® 
a x Px ’ (57) 
x 


suggesting a next approximation on from Eq. (52) is 


F 
ows | yptsin ae (58) 
Ox 


Eq. (58) follows from the integration of Eq. (56). This pro- 
cedure can be iterated to an n’” order WKB approximation 
in the boson FW interpretation. 


9. 2nd order relativistic WKB approximation 
Expanding ® (x) in powers of h 


@ (x) = By (x) +hH®, (x) + H°D, (x) +... (59) 
then substituting Eq. (59) into Eq. (56) gives 

(@/, + no! +...)” — ih (OY +n" +...) = p2. (60) 
Equating terms in powers of h, in Eq. (60), I see 


+O, =p,, 20/0! —io! =0, (61) 


so that 


® = +f rax, ®, = 5 log, ®,. (62) 


Eqs. (62), Eq. (59), and Eq. (53) give 


ey (63) 


This result is identical to the Klein-Gordon WKB approxi- 
mation because on the way I approximated FW Eq. (18) with 
Eq. (16). 


10. Discussion of fermion and boson WKB 
approximations 


When I made my WKB approximation in Section 7, from 
what I will show in Section 14 and Appendix B to be the 
correct form Eq. (18) of the FW equation, I saw in the non- 
exponential factor of Eq. (48) appearances of p; changed to 
appearances of a; going from boson to fermion descriptions. 
Ihave noticed, that for tunneling through slowly varying po- 
tentials, Eq. (48), Eq. (50) and Eq. (63) give the implication: 


Eq. (64) suggests that the WKB approximation for fermions, 
in the correct FW representation, will take the form 


Dg (65) 


Going from Eq. (63) to Eq. (65) I swapped an appearance of 


p; to +a, from boson to fermion descriptions, as suggested 
by Eq. (48), Eq. (50) and Eq. (63). 


11. Connection formulae 
It is suggested by Section 8 that for the WKB approxi- 
mation to be appropriate 


Op, 
Ox 


Ip, | > lh (66) 


Eq. (66) suggests that when p, ~ 0 the WKB method fails. I 
need a connection formula between the regions where ae # 
0. 

I assume that when Eq. (66) fails 


me > |px| . (67) 


so that I may treat the connecting region with the c = oo 
non-relativistic limit of the FW Eq. (18). The wavefunction 
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in the connecting region is then given by the Airy function 
solutions to the Schrédinger equation. This approach leads 
to the two connected solutions for ¥ between the classical 
and non-classical regions: 

Classical region p” > 0 and x < x, 


i i 
exp (+5 J, |p| dx’) exp (-5 Ms Ip| dx’) 
A B 
v IPI 


Non-classical region p” < 0 and x > x, 


1 
exp (+5 es |p| dx’) 
+B, , (69) 


viol 


2A, = -i (A, +iB,) exp ( +i2/4) (70) 
and 
By = +i (A; —iB,) exp ( - iz/4). (71) 


In Eq. (68) and Eq. (69) I have used the asymptotic forms of 
the Airy functions and taken p = 0 to be at x,. 

The wavefunctions Eq. (68) and Eq. (69) are in the WKB 
form derived for relativistic fermions in Section 7 to Section 
10, therefore the solutions in both regions can be extended 
to relativistic energies. 


12. Extention to periodic structures 


Consider the periodic potential 


iGx ) (72) 


V(x) = Vo Veexp (+= 
G 


Eq. (72) is a Fourier expansion of the potential, where G = 
2nn/a and n is positive and negative integers. The wave- 
function follows from Bloch’s theorem for periodic crystals 


GO = exp (+=) > Ag exp (-#*), (73) 
& 


where g = 2am/a and m is positive and negative integers. 
Substituting Eq. (72) and Eq. (73) into 


[FL - (E+me*+V)| pH = 0 (74) 


and multiplying through by exp (i(G’ — k)x/h), where G’ = 
2an’ /a, then integrating x over one period a, I find that the 
first, second, and third terms are only non-zero for g = G’ 
while the fourth term is only non-zero for g—G = G’. There- 
fore, with rearrangement 


DY, (6g6"W — Vegi) Ag = (Ex me?) Agr, — (75) 
& 


where 


W, =- CS a at?) (ck — cg) . (76) 
The eigenvalue problem described by Eq. (75) can be 
solved for varying k to give the eigenvalues (E + me’). The 
relationship between k and E defines the band structure. 
For planar structures of period D I have by Bloch’s the- 
orem 


P(x + D) = exp (“2 )H a, (77) 


which in terms of transfer matrices can be written as 


M¥ = exp (“> )# = a9 —> |M —Al| =0, (78) 
where 
M M 
M= 11 ) . 719 
Ge Mn (7) 


For fermions to propagate through the periodic structure K 
must be real, therefore by Eqs. (78) 


M,,+M. 
-1 <cos(*2) = no <4. (80) 


Band-gaps arise in 1D periodic potentials if the condition in 
Eq. (80) is not met. 


13. Calculation of transfer matrix parameters 


13.1. Calculation of the angle of refraction 
To match the Dirac spinors of Eq. (8) and Eq. (9) either 
side of a potential step at x = O it must be that 


P| -0 = Plo» (81) 


therefore, since 
p=xXpcosd+t ppsing, (82) 
l arrive at Snell’s law for the relativistic fermion 


sin pO _ p® 
sing pO” 


(83) 


I do not analytically verify non-normal incident relativistic 
tunneling and transmission in this article. 


13.2. Tunneling between layers of varying electron 
mass 
In free space by special relativity 


Bl = me! (1 rn (‘)) (84) 


inside a layer by special relativity 


(E+V) =m"?cr4 (1 +(4 )) , (85) 


c* 
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where m* is the effective mass of the electron in the layer and 
c* is the corresponding effective speed of light. Comparing 


Eq. (84) and Eq. (85), I find 
me? 2g(k/c)V +( V 


z=(8) + 
k2 \ m*c* m*ct2— m*k2 


(me? — m*c*?) (me? + m*c*?) 1 
$+  — —. (86 
(m*c*)? ke? on 


In the last term of Eq. (86) I used the difference of two squares 
to improve numerical accuracy. When |k| < c the function 
g(k/c) is given by the Taylor expansion 


k 2 


2 4 
g(k/c) = 1+(<) aggro _ (k/c) 


2 8 


+... (87) 


and when |k| > c the function g(k/c) is given by the Taylor 


expansion 
ky? ky 
g(k/e) = (*) y 4 CD" Cs). ee 
c 2 8 
Computationally I add smallest terms together first for im- 
proved numerical accuracy. In terms of g(k/c) and the x- 


component of electron momentum m*q, I calculate the pa- 
(+) 


rameter a,” to be 
29 (k V +tmi*c*? 
a = —me g(k/c)+V+m*c . (89) 
q cm" ¢, 


13.3. Tunneling into a layer where electron mass 
vanishes 
In free space the fermion momentum dispersion is given 
by Eq. (84) whereas inside the layer electron mass has van- 
ished and therefore 


E+V =qce", (90) 


where q is the fermion momentum in the layer, and c* is the 
speed of light in the layer. Comparing Eq. (84) and Eq. (90), 
I find 


qc* = me7g(k/c)+V (91) 
and 
a,=-~. (92) 
qx 


14. Results 


This section is restricted to quantum tunneling. For the 
structures I investigate in this section I introduce a length 
scale a given in terms of electron mass m, the speed of light 
c,andh 

h 
a=5—. (93) 


mc 


Tansmission T 


Abs. error in T 


: p 
Incident mec 


Figure 1: (Color online) (a) Transmission through the poten- 
tial barrier Eq. (94) as a function of incident momentum p, for 
fermions (red) and bosons (orange). (b) Convergence of the 
fermion transmission with increasing number of relativistic cor- 
rection Taylor expansion terms 0, 5, 10, 20, 40, 80 as labelled. 


In Appendix B I investigate reflection of electrons from a 
sharp step in potential. In Appendix B I find exact analytic 
agreement for R between the FW and Dirac spinor represen- 
tations. 

The striking feature of all my results is that I see exact 
analytic agreement between T and R calculated in the FW 
and Dirac spinor representations for all momentum p, even 
when |p| > mc. This is a significant finding since when 
|p| > mc the Taylor expansion of £ is a divergent series. I 
think that conservation of energy is restricting my evaluation 
of £ on ¥©) to the particular values required to meet the 
physical boundary conditions for the Dirac Hamiltonian. 

I model a resonant-tunneling diode to compare my WKB 
approximations of Section 9 and Section 10 as the final ex- 
ample in this section. 


M. B. Doost: Preprint submitted to Elsevier 


Page 7 of 14 


Foldy-Wouthuysen WKB method for Dirac tunneling 


Tansmission T 


— Dirac 
—— Klein-Gordon 


Abs. error in T 


p 


Incident me 


Figure 2: (Color online) Repeat of Fig. 1 except that now elec- 
tron mass vanishes inside the barrier. 


For my numerical computations I use Python 3.8. 


14.1. Tunneling through a barrier 
The system under consideration is a potential barrier de- 
scribed by 


vay=-f ae 


3mc 


for |x| >a, 


for |x| <a. OP 


Since V is negative, the momentum of the electron inside the 
barrier becomes imaginary when the following condition is 
met: 


—mc* <E+V <+mc’, (95) 


then the electron is quantum tunneling. 
The analytic transmission for Eq. (94) is derived in Ap- 


pendix D as: 
4K? 


its — — _ Fe 
") 2... [ 2mqa 
4x? + (1 — «?)* sin 5 


(96) 


where I take « to be x, in the FW representation, and «, in 

the Dirac spinor representation. For the derivation of Ap- 

pendix D in the FW representation 
(+) 

a) % — 4 E+ me? 


kK = = 97 
f a kK, E+mce2?+V oa 


It was shown by the authors of Refs. [14, 15] by deriving T 
in the Dirac spinors representation of quantum tunneling 


(4) Ge Et me 
K = 


d k, Etme2+V’ 8) 
which is in exact analytic agreement with oes as expected. 

The numerical results demonstrating relativistic tunnel- 
ing of electrons are shown in Fig. 1. For comparison I also 
show the transmission of tunneling bosons with identical 
mass to the electron. 

Fig. 1(b) shows the absolute error calculated as the dif- 
ference between T evaluated using the Taylor expansions 
g(k/c) in x or the Python 3.8 cmath module square root 
function. 

The region of slow convergence around |k/c| = 1 in 
Fig. 1(b) is an artifact of the point in g(k/c) about which I 
have taken my Taylor expansions. 


14.2. Tunneling through a barrier with vanishing 
electron mass inside the barrier 
In this subsection I almost repeat the demonstration of 
Section 14.1 except that now inside the barrier the electron 
mass is strictly vanishing. 
Since m = 0 when |x| < aI must re-evaluate x, for 
Eq. (96): 


2 
ges 4x E + mc (99) 
d Dy E+V 


which is in exact analytic agreement with c®, as expected. 

The numerical results demonstrating tunneling of rela- 
tivistic electrons are shown in Fig. 2. For comparison I also 
show the transmission of tunneling bosons with identical 
mass to the electron. 

In Fig. 2(b) the absolute error is calculated as the dif- 
ference between T evaluated using the Taylor expansions 
g(k/c) in x or the Python 3.8 cmath module square root 
function. 

The region of slow convergence around |k/c| = 1 in 
Fig. 2(b) is an artifact of the point in g(k/c) about which I 
have taken my Taylor expansions. 
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Potential V/mc? 


Transmission Tg 


—0.02 


Applied potential AV/mc? 


Figure 3: (color online) (a) Potential profile for a representa- 
tion of a resonant-tunneling diode under a range of bias po- 
tentials AV /mc? 0.00, 0.25, 0.50, 0.75, 1.00 as labelled. (b) 
Transmission as a function of bias potential AV for a series of 
left hand side electron incident energies E/mc? 0.400, 0.425, 
0.450, 0.475 as labelled, with I, taken to be +a. (c) Dif 
ference in transmission, as a function of bias potential AV, 
between using +a‘ or p, for T;,. 


i 


14.3. Transmission through a resonant-tunneling 
diode using the WKB approximation 

I calculate the transmission through a resonant diode un- 
der a range of bias potentials AV. I show a series of poten- 
tials describing increasing AV in Fig. 3(a). The transmission 
as a function of AV is given in Fig. 3(b) for a series of inci- 
dent energies. I calculate transmission with the two different 
WKB approximations, derived in Section 9 and Section 10, 
for comparison in Fig. 3(c). 


The resonant-tunneling diode that I consider includes a 
planar cavity with a smoothly varying potential. For the 
wavefunction to traverse the cavity the employment of WKB 
transfer matrices developed in Section 9 and Section 10 are 
required: 


exp (-2 pdx) 

— 0 
V ofr, . = 

exp (+5 7 pdx) 


b 


Ql 
iS) 
So 


0 


(100) 


I; is p; or +a depending on whether I use FW Eq. (16) or 
Eq. (18) inside the cavity. The subscript 0 denotes evaluation 
at xy = —a/2 the left hand side of the smoothly varying 
region for which the transfer matrix is applied, | is atx; = 
+a/2 the right hand side of the smoothly varying region. I 
define Gg, @;, bo, and b, in the same way as Section 5. 

Since I have established in Section 14.1 and Section 14.2 
that the Taylor series for g(k/c) gives convergence to the 
square root function of the Python 3.8 cmath module, I dis- 
pense with the Taylor series and use the cmath square root 
function to calculate the results in this subsection. 

I see in Fig. 3(b) and (c) that as the energy of the inci- 
dent electrons moves to higher values the first resonance of 
the diode moves towards lower values of AV, as expected. 
I observe in Fig. 3(b) negative transmission as expected ac- 
cording to the Klein paradox [14, 16]. 

In Fig. 3(b) and (c) I see no qualitative difference be- 
tween my two WKB approximations, showing that exponen- 
tial decay is the dominant factor in the cavity for this exam- 
ple. 


15. Conclusion 


In this work I have developed Dirac fermion transfer ma- 
trix methods, for all energies, in the FW representations. 
I have also developed WKB approximations, to all orders, 
and discussed the validity of the approximations. I have in- 
troduced connection formulae between regions of real and 
imaginary fermion momentum. 

I have shown the extension of my method to 2D periodic 
structures for band engineering [1, 2, 3]. 

I have alerted the reader to the limits of my methods by 
revealing a GF paradox that occurs when transforming be- 
tween Dirac and FW representations. 

I have verified my methods for massive and effectively 
massless fermions of all momentum by analytic compari- 
son with the Dirac spinor representation of fermion tunnel- 
ing. I have shown that the FW representation is an exact 
description of Dirac fermions and is not restricted to semi- 
relativistic energies as first envisioned [11]. 
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A. Discussion of commutation 


In this appendix I will provide a discussion of 
L=-Lep, =—cp,L. 


Consider the form of the operator £ 


a 72 a 74 ~ 716 
L=mce? fe (ie pe + a 
2 |mc 8 | mc 16 | mc call: 


(102) 
where 
2 

1| p hn? |[d* ad 
cM peat | ge | 103 
2 =| 2m2c2 FE 7 dy* a 
ifpl_ nt [a d‘ d‘ 
= #9 z | . a0 
8 | mc 8m4c4 | dx4 dy*dx2  dy4 


and the series of polynomial differential operators continues 
up to n'” order. It is clear that I can write 


qutm d qutm-1 qutm-1 d 
pee ne 105 
dx"dym dx dx"-lqym dx"-lqym dx ( ) 
However it is also the case that 
n qd’! n—1 d 
— fay) = = ed aes - cee (106) 
dy" dxdydy"-1_— dy"! dy dx 


when every term in f is a function of not just y but also x. 
I see from Eq. (105) and Eq. (106) that I can extract d/dx 
from either side of the operator £, in these cases. 

When f(x, y) = f(Q), then in general 


d"” 
dy" 


fM #09, (107) 


but 
d 
ae (108) 
x 
therefore in general 
d" d"! dx d 
0 SS =0. 109 
Pap OF aia (109) 


In the case described by Eq. (109) I cannot extract d/dx 
from either side of the operator d"/dy". 

Eq. (101) is true when I assume that £ acts on a function 
of not just y but also x. 


B. Reflection at a step potential in 
Foldy-Wouthuysen and Dirac spinor 
representations for comparison 


Consider a fermion incident on a step at x = 0. The 
wavefunction for x < 0 is 


2 ip,x 2 ip,x 
HE exp (+ Px +B fame exp [ — Px ; 
—cp, h +cp, h 


(110) 
and for x > 0 
E+V+mc* iq, 

F( ey ) exp (+ i ) ; (111) 
By continuity of wavefunction at x = 0 I have 

(1 - B)cp, = Fe*q,, (112) 

(1+ B)(E+me*) =F(E+V+me"). (113) 
Dividing Eq. (112) by Eq. (113) gives 

= 2 
BSF ee Ee ge). (114) 
1+B pyc E+Vtme 4 


Now I turn to the FW representation. The corresponding 
wavefunction for x < 0 is 


exp (+2) + Bexp (-2=), (115) 
and for x > 0 
Fexp (+4). (116) 


Applying my boundary condition at x = 0, continuity of 
wavefunction and continuity of wavefunction with respect 
to £, I obtain 


(1- B)a® = Fa, (117) 
1+B=F. (118) 
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Dividing Eq. (117) by Eq. (118) gives 


a A 
1-B 1 Px¢ E+V+4mc* (+) 
1+B yg  q.c* E+mce? f 
0 
In Section 3 I derived R = | B|?, therefore 
taer 
= ; 120 
| l+«K ge) 


where I use —_ for the Dirac equation and a for the FW 


equation. Substituting either x or cS into Eq. (120) gives 


analytically identical T and R between the Dirac and FW 
representation, as expected. I also calculated oa in the case 
of m vanishing for x > 0 to find exact analytic agreement 
in T and R between the Dirac and FW representations, as 
expected. 

x corresponds to a), and is the mathematical conse- 
quence of there being two forms for the FW wavefunction 
Eq. (18). The FW wavefunction describes the electron po 
or positron ¥) [11]. 

IT also calculate, for useful comparison, the transmission 
of bosons with the following Klein-Gordon result 


_ 4,c* 


~ pe 


Kkg (121) 


C. Examination of a massless fermion 
tunneling through a magnetic delta barrier 
in the Foldy-Wouthuysen representation 
Consider a massless fermion tunneling through a mag- 

netic delta barrier B = Bod (x) 2, where B = V X A. The 

barrier is described by A = 0 except for A = A,j # 0 when 


x > 0. The Dirac spinor wavefunction is a solution of 
[ca -(p-—A)]V= EW. (122) 


Using the approach of Section 2, Eq. (122) can be trans- 
formed to 


lev (p- Ay ye = pp | (123) 
Taylor expanding Eq. (123) around p = 0 gives 
[FM - (E+cA,)| P =0, (124) 


where M is the corresponding linear differential operator. 
Using the approach of Section 4 I calculate, that at the bar- 
rier, the wavefunction is continuous and also continuous with 
respect to M, where I define 


ExtcaA,~ 


MEH = + pe) | (125) 


CPx 


I now return to the Dirac spinor representation. In the 
region A = A,j # 0, the general solution for p = +q,X is 
given by [22, 23, 24] 


+E igyX 
= ; exp | + ; 
a 7 ) ( h ) 


(126) 


where cq, = ,/E* —c? A}. The reflection is calculated by 


matching wavefunctions at the step in vector potential, sim- 
ilar to Appendix C for a step in electrostatic potential. 

Making use of my boundary conditions I find the reflec- 
tion coefficient R takes the form of Eq. (120), with 


* q,—iA 
Ky= cai g (127) 
Cc Pp, 
for the Dirac spinor representation and 
Cc Et c*A 
geo (128) 
c*dy E 


for the FW representation. p = p,X = (E/c)X when A = 0. 

The conflict between Eq. (127) and Eq. (128) shows that 
R calculated in the FW representation does not reproduce R 
calculated in the Dirac spinor representation, FW wavefunc- 
tions cannot be used for the calculation of tunneling through 
magnetic barriers. 


D. Analytic Green’s functions in 1-dimension 


Consider 
[FL - (V+ E+mce’)| G® (x, x’) =5(x—x’) . (129) 


Let ¥ y and v7 r be the solutions of Eq. (18) which separately 
satisfy the boundary conditions at x = —a and x = +a. The 
GF of Eq. (129) for the case of a 1D barrier can be written 
as 


Pe (x) PS (x) 


(+) 1) 
G (ox ) = Vo : (130) 
where x. = min (x, x’) and x, = max (x, x’) and 
WH = FPP icnlLP® + PD ich ly (131) 


is the Wronskian, which does not depend on x. 
To prove Eq. (130) and Eq. (131), note that by Eq. (129) 


x! + 


€ 
+ lim LG (x)dx=1, (132) 
e—>' / 


XxX, SE 


which implies 
+ichLG™ (x,,x') + ichLG™ (x_,x’) =1, (133) 


then, applying £ to Eq. (130) written out explicitly 


= rn _ 1 LY, (x) P p(x’) for x <x’, 
ee ek { P, (x! LY¥R x) for x>x’, 
(134) 


and substituting into Eq. (133) gives Eq. (131). 
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I now calculate the GF for a barrier of width 2a. 


egal imkx 
h 
im* qx im* qx 
Cexp {+ + Dexp | -——— |x| <a, 
h h 
Aexp (+ an + Bexp( - mk x<-—a. 


(135) 


x > +a, 


Continuity of ‘¥; (x) and continuity of ‘Y; (x) with respect to 
L£ gives 


ag + ay 


imka im*qa 
— , (136 
2a, exp (+ h ) exp ( A ) (136) 


a, — a ; im* 
D=— exp (+ tt) exp (+ #). as 


C= 


2a, 


¥, the left hand solution is given by (+x) a W,(-x). 
Therefore the Wronskian is given by 


W = 2icha, (C? — D?) . (138) 
I evaluate Eq. (130) to give G(+a, —a): 


IQ 


am* am* . 
ich (a? + a?) sin ( “) — 2cha,a, cos ( “a 


(139) 


The transmission for a barrier embedded in a vacuum 
is given as the ratio of G(+a,—a) and the free space GF 
GfS(+a, —a): 


2 


BGH 2 |. a6) 


2: agate 
T(k) = |2cha,G (+a, —a) | Fas (aa 


The free space GF can be written as 
[zL - (E +mc’)| G® (x, x’) =6(x-x’) , (141) 


therefore by symmetry, when x’ = 0, 


(-“*) for x <0, 


(4) _ lt 
G(x, 0) = A* ke 
xp (+ i ) for x >0. 


(142) 


Integrating Eq. (141) across [O_, 0 Al gives 


= imkO - imkO_ 
FichA™L exp (+ 7 us )iena2 exp (-“*) =] 


(143) 
and letting 0, > O and 0_ — 0,1 find 
i: ee (144) 
2icha*? 


E. Green’s function paradox 


I reveal a GF paradox occurring when transforming be- 
tween the FW and Dirac representations. 


Let 
[H -(E+V)G(rr’) =6(r—-r') 1, (145) 
therefore 
eS LH -(E+V)Je'Set'§ [¢(r,r’)| e'® 
= et fe'S§(r—r’'). (146) 
However since 
eS fe 'S5 (r-r’) = 16 (r—-r’) , (147) 


and since I have already seen in Section 2 


etis [H -(E+ V)\e75 _ [FL - (V + E+mc’)| i 
(148) 
Ihave 
8 [G (rae = G(r) 
= eG (rr) et =G (rr) 1. (149) 


Eqs. (149) imply that the valid GFs of the Dirac equation are 
diagonal. Now consider the identity 


[H —(E+V)][H +(E+V)|K (r,r’) = 15 (r-7') , 


(150) 
so that 
[H+(E+V)|K (7,1) =G(r,7) . (151) 
Multiplying out the brackets in Eq. (150) I obtain 
242102 p\? / ! 
-'n |w?+ (2) 1K (rr) =8(r—r'). (152) 
I see from Eq. (152) that in 1D free space (V = 0): 
ip |x — x’ 
cn? K (x,x") = ee (153) 
? 2ip ? 
in 2D free space: 
a2  — ia (Ple-e'| 
h°K (p,p) = ~H ———— }s 154 
(p, 0") re ( (154) 
and in 3D free space: 
(7 |r—r'| ) 
os ores 
272 / 
h’K (r,r ) = —————_—_. 155 
WK (1) = (155) 
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Substituting the identity Eq. (151) into Eq. (149) 
eS (H+ (E+V)e SK (r,r') = G™ (r,r’) 1 (156) 
and so 
| VirCP + mck + (E+ v)| K (r, a = G® (r, r’) i. 
(157) 


By inspection of Eq. (157) and the energy momentum rela- 


tion (E+ V) = V/p2c? + mc4, 


G® (rr!) =2(E+V)K (r,7’) (158) 
and substituting Eq. (158) into Eq. (149) I also have 
G(r,r’) =2(E+V)K(r,r’). (159) 


Eq. (159) is consistent with the relativistic correction to the 
1st Born approximation in 3D. 

Let us test Eq. (158) and Eq. (159) by making an analytic 
comparison between the GF for 1D free space Eq. (142) and 
the GF composed of Eq. (159) and Eq. (153): 


and so I have an apparent paradox. 
Let us attempt to resolve the paradox. I have shown in 
Appendix D that there are two possible forms of the FW GF: 


(161) 


for an electron delta source (—) or a positron source (+). If I 
add my two expressions G)(x, x’) I find 


GO ( =) +G® (3 x) =2EK ( x’) : (162) 
However if I try to prove that in general 
G(r,r') = GO (7,7) +6 (717), (163) 


I arrive at a contradiction. The Klein-Gordon GF has the 
continuity conditions of the Klein-Gordon equation, (E+V) 
can be a discontinuous function, therefore 2(/E +V)K can be 
discontinuous at boundaries. Therefore Eq. (159) implies G 
may be discontinuous. G has the same boundary conditions 
as , Eq. (159) may violate the boundary conditions for G. 

Let us now discuss a possible cause of the paradox. To 
evaluate the FW transformations defined to be 


et!S — cosO + Ba- 


|» 


sin 6, (164) 


» 


I am required to know the momentum / also at the delta 
source. To make this calculation of f at the delta source 
consider Eq. (142) from Appendix D, 


imkx 


1 exp (- ) for x <0, 


exp h 


G(x, 0) = (165) 


) for x >0. 


I see from Eq. (165) that 
Dx (x = 0) G(x, 0) = Jim, pG (0,.,0) = +mkG (0,0) . 
(166) 


Eq. (166) indicates that when x = x’, the momentum is 
undefined. Therefore at the delta source the FW transforma- 
tion is undefined. In 1D it might be impossible to evaluate 
e*!S at the delta source. 
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